clear; clc; close all;

%% 参数设置
k = 3;
mu = pi;                 % 固定 mu = pi
h = 0.01;
T = 600;
transient = 550;
N = T / h;
N_trans = transient / h;

eta_list = linspace(-3*pi, 3*pi, 200); % eta 变化范围
num_LE = 5;
LEs_all = zeros(length(eta_list), num_LE);

%% 并行计算李雅普诺夫指数
parfor i = 1:length(eta_list)
    eta = eta_list(i);
    x0 = [1e-9; 0; 0; mu; eta];  % 初始状态，mu固定，eta变化
    LEs = computeLEs(k, h, N, N_trans, x0);
    LEs_all(i,:) = LEs(1:num_LE);
end

%% 绘图
figure; hold on; box on;
plot(eta_list, LEs_all(:,1), 'r-', 'LineWidth', 1);    % 最大LE
plot(eta_list, LEs_all(:,3), 'g-', 'LineWidth', 1);    % 最小LE
plot(eta_list, LEs_all(:,4), 'b-', 'LineWidth', 1);    % 最小LE
plot(eta_list, zeros(size(eta_list)), 'k--', 'LineWidth', 1);

xticks(-3*pi : pi : 3*pi);
xticklabels({'-3\pi', '-2\pi', '-\pi', '0', '\pi', '2\pi', '3\pi'});

xlabel('\mu','FontSize', 12);
ylabel('LEs', 'FontSize', 12);

legend({'LE_1', 'LE_3', 'LE_4'}, 'Location', 'best');
grid on;

%% 你之前的 computeLEs, system_ode, jacobian_matrix, RK4_step 函数保持不变


%% 子函数：计算李雅普诺夫指数
function LEs = computeLEs(k, h, N, N_trans, x0)
    % 系统维数
    n = length(x0);
    
    % 初始化轨道和扰动向量
    X = x0;
    Q = eye(n); % 正交扰动向量矩阵
    
    % 李雅普诺夫指数初始化
    sum_LE = zeros(n,1);
    
    dt = h;
    steps = N;
    
    for step = 1:steps
        % 主轨道积分一步（RK4）
        X = RK4_step(@(x) system_ode(x, k), X, dt);
        
        % 线性化方程积分扰动向量（雅可比矩阵）
        J = jacobian_matrix(X, k);
        Q = Q + dt * J * Q; % 欧拉法积分扰动，简单可改为RK4
        
        % 正交化 QR 分解
        [Q, R] = qr(Q);
        
        % 李雅普诺夫指数累积（对角线）
        log_diag = log(abs(diag(R)));
        
        if step > N_trans % 丢弃暂态
            sum_LE = sum_LE + log_diag;
        end
    end
    
    % 计算平均李雅普诺夫指数
    LEs = sum_LE / (dt*(steps - N_trans));
    LEs = sort(LEs, 'descend')';
end

%% 子函数：系统方程
function dx = system_ode(x, k)
    % 状态变量 x = [x; y; z; u; v]
    dx = zeros(5,1);

    dx(1) = x(2) + x(3) - k * cos(x(5)) * x(2); % \dot{x}
    dx(2) = -x(1) + x(3);                        % \dot{y}
    dx(3) = -x(1) - x(3) + k * cos(x(4)) * x(1); % \dot{z}
    dx(4) = x(1);                               % \dot{u}
    dx(5) = x(2);                               % \dot{v}
end

%% 子函数：计算雅可比矩阵
function J = jacobian_matrix(x, k)
    % 这里的 x(4)=eta，x(5)=mu，注意论文中对应位置含义
    eta = x(5); % 论文中 eta 和 mu 变量位置可能和代码中对应不同，确认清楚！
    mu = x(4);
    
    J = zeros(5,5);
    J(1,2) = 1 - k * cos(eta);
    J(1,3) = 1;
    
    J(2,1) = -1;
    J(2,3) = 1;
    
    J(3,1) = -1 + k * cos(mu);
    J(3,3) = -1;
    
    J(4,1) = 1;
    J(5,2) = 1;
end
%% 子函数：RK4一步积分
function x_next = RK4_step(f, x, h)
    k1 = f(x);
    k2 = f(x + h/2 * k1);
    k3 = f(x + h/2 * k2);
    k4 = f(x + h * k3);
    x_next = x + h/6*(k1 + 2*k2 + 2*k3 + k4);
end


